function ir = image_resize(im, h, w)
% image_resize  - resize a HxWxD numeric image to new measures
%
% FORMAT:       ir = image_resize(im, h [, w])
%
% Input fields:
%
%       im          HxWxD image (height x width x color depth)
%       h           new height (or maximum dimension)
%       w           new width
%
% Output fields:
%
%       ir          resized image
%
% Note: if both h and w are given, one of the two can be set
%       to 0 which will lead to auto-detection of the second
%       parameter; alternatively only h can be given which will
%       set the larger of the two dimensions to this size

% Version:  v0.7g
% Build:    9100717
% Date:     Oct-07 2009, 5:00 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   ~isnumeric(im) || ...
    isempty(im) || ...
    numel(h) ~= 1 || ...
   ~isa(h, 'double') || ...
    isinf(h) || ...
    isnan(h) || ...
    h < 0 || ...
    h ~= fix(h) || ...
   (nargin < 3 && ...
    h < 4)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
imc = class(im);
if nargin < 3
    maxsize = true;
    w = 0;
else
    maxsize = false;
end
if numel(w) ~= 1 || ...
   ~isa(w, 'double') || ...
    isinf(w) || ...
    isnan(w) || ...
    w < 0 || ...
    w ~= fix(w) || ...
    all([h, w] < 4)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad argument.' ...
    );
end

% get image dims
ims = size(im);
if numel(ims) < 3
    ims(3) = 1;
end

% get requested size
if maxsize
    [ms, msp] = max(ims(1:2));
    if msp == 1
        newheight = h;
        newwidth = round(h * ims(2) / ms);
    else
        newheight = round(h * ims(1) / ms);
        newwidth = h;
    end
else
    newheight = h;
    newwidth = w;
    if h == 0
        newheight = round(w * ims(1) / ims(2));
    elseif w == 0
        newwidth = round(h * ims(2) / ims(1));
    end
end

% prepare output
ir = zeros(newheight, newwidth, ims(3));

% get interpolation points
imht = ims(1) / newheight;
imhs = 0.5 * (imht + 1);
imhe = (ims(1) + 1) - imhs;
imwt = ims(2) / newwidth;
imws = 0.5 * (imwt + 1);
imwe = (ims(2) + 1) - imws;

% smooth first?
if any([imht, imwt] > 1)
    if ims(3) == 1
        im = im(:, :, ones(1, 3));
    end
    im = smoothdata3(double(im), [imht, imwt, 0], 0.01);
    if ims(3) == 1
        im = im(:, :, 2);
    end
end

% interpolate each color dim
for c = 1:ims(3)
    ir(:, :, c) = flexinterpn_method(im(:, :, c), ...
        [Inf, Inf; imhs, imws; imht, imwt; imhe, imwe], 0, 'cubic');
end

% make class correct
switch (lower(imc))
    case {'int16'}
        ir = int16(fix(max(min(ir, 32767), -32768)));
    case {'int32'}
        ir = int32(fix(max(min(ir, 2^31 - 1), -(2^31))));
    case {'int8'}
        ir = int8(fix(max(min(ir, 127), -128)));
    case {'single'}
        ir = single(ir);
    case {'uint16'}
        ir = uint16(fix(max(min(ir, 65535), 0)));
    case {'uint32'}
        ir = uint32(fix(max(min(ir, 2^32 - 1), 0)));
    case {'uint8'}
        ir = uint8(fix(max(min(ir, 255), 0)));
end
